/**
 * @file   main7.cpp
 * @author Carlos Hernández García <carloshergar@usal.es>, Camilo Ruiz, Enrique Conejero
 * @date   Fri Apr  10 2009  
 *
 * Copyright (c) 2009 Carlos Hernández García
 *
 * @brief  Main program, in which we have running correctly the routines Charge_deposition and Current_deposition
 *
 * 
 * Carlos Hernández García                          \n
 * <carloshergar@usal.es>                           \n
 * Universidad de Salamanca                         \n
 * Plaza de la Merced s/n                           \n
 * E-37008 Salamanca, SPAIN                         \n
 */
/********************************************************************************/


#include <iostream>
#include <vector>
#include <string>
//#include <cstring>

#include <alloc.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <complex>
#include <fstream.h>

#include "../include/Parameter.h"
#include "../include/carpetas.h"
#include "../include/gridobject.h"
//#include "../include/laserfield.h"
#include "../include/routinesO.h"
#include "../include/constants.h"


#define complex complex<double>

using namespace std;
using std::cout;
using std::vector;
using std::string;

Parameter<string> jobname("jobname","Run1");
Parameter<double> lambda_nm("lambda_nm",32, "Wavelength in nm");
//                                          ^^ optional comment
Parameter<vector<int> > point("point", vector<int>(3,0) );

Parameter<int>    np0("np0",4000,"Particle number"); 
Parameter<int>    ndg0("ndg0",8192,"Spatial grid points"); 
Parameter<int>    ndt0("ndt0",32768,"Temporal grid points");

//Mode of calculation
Parameter<int>    FieldsCalculation("FieldsCalculation",1,"Mode of calculation of the fields (1-Maxwell, 0-Potentials)");
Parameter<int>    ChargeCalculation("ChargeCalculation",0,"Mode of deposition of the charge (1-Lichters, 0-Enrique)");
Parameter<int>    CurrentCalculation("CurrentCalculation",0,"Mode of deposition of the current (1-Lichters, 0-Enrique)");

//Field parameters
Parameter<double> w0("w0",0.05,"Laser frequency");
Parameter<double> E0("E0",4,"Laser Amplitude");
Parameter<int>    nw("nw",20,"Number of laser cycles");

//Density
Parameter<double> density("density",15.,"Density in critical density units");

//Parameter for the average
Parameter<int> na("na",3,"Smoother");

//Output folders
//Parameter<string> foldername("foldername","pic");
Parameter<int> dir_limit("dir_limit",50,"Maximum number of directories that can be created");





int main(int argc, char *argv[])
{

  // read configuration file
  ConfigFileParser cfg( argv[1] );
    
  ParameterMap pMap;
  // initialize all Parameter<T> type objects from the config file
  pMap.init( cfg );
  
  
  cout<<"Wavenumber: "<< 1/( lambda_nm * 1e-9 ) <<"\n";
  
  // string/map access:
  double x; cfg.getValueInto(x, "lambda_nm"); // ConfigFileParser style
  double y = pMap["lambda_nm"];               // ParameterMap style
  string s = pMap["jobname"];                 // kind of "chameleon"
  
  try {
    cout << "Trying to pass a Parameter<double> to an int:\n";
    int i = pMap["lambda_nm"];
  } catch ( std::runtime_error e ) {
    cerr << "   ->    " << e.what() <<"\n";
  }
  
  cout<<pMap;
  
  
  
  /*************************************************/
  //Create folders and files.
  FILE *fp0; //Fichero de prueba
  FILE *fo1, *fo2, *fo3; //Ficheros de salida de datos
  //string hh=pMap["foldername"];
  //char nomcarpeta=hh.c_str();
  
  //Creating a new data directory
  
  if(CrearCarpeta(dir_limit)==0)
    {	
      printf("Has superado el limite de directorios que tu pusiste, cazurro!!");
      exit(1);
    }
  
  
  fo1=fopen("output1.txt", "wt+"); //Fichero donde se almacenan los campos.
  if(fo1==NULL)
    {
      printf("Error opening the file.\n");
      exit(2);
    }
  fo2=fopen("output2.txt", "wt+"); //Fichero donde almaceno las cargas.
  if(fo2==NULL)
    {
      printf("Error opening the file.\n");
      exit(2);
    }
	fo3=fopen("output3.txt", "wt+"); //Fichero donde almaceno las cargas.
	if(fo3==NULL)
    {
      printf("Error opening the file.\n");
      exit(2);
    }
  
  
  /*************************************************/
  // Define variables
  /*************************************************/
  
  
  
	//int iaux,npion,npim;
	int npion=0;
  
	int np=pMap["np0"];//4000; //40000;//400;//3; //4000;              ///< Particle number
	int ndg=pMap["ndg0"];//1024*8; //1024*8; //200; //1024*8;           ///< Number of points in the spatial grid
	int ndt=pMap["ndt0"];//1024*32; //1024*32; //5000;//1024*32;      ///< Number of points in the temporal grid
  
  
	//Variables for output
	int nscleft, nscright;
	//int if1,if2,if3,if4,if5;
	int nsc;
	int idj[10000];
	
	
	//Field parameters
	//double w0=0.05; //0.057;
	//int nw=20;
	//double E0=4.;
	laser laser0(pMap["w0"],pMap["nw"],pMap["E0"]);
	
	//Calculation parameters.
	double lg=(laser0.nw+1)*laser0.lambda/4; //WATCH Fix here to make deltag=c*deltat
	double tmax=(laser0.nw+1)*2*Pi/laser0.w0;
	double deltag=lg/ndg;
	double deltat=tmax/ndt;
	
	//Material parameters
	double thickness=laser0.lambda;
	double xgleft=lg/2.-thickness/2.;
	double xgright=xgleft+thickness/2.;
	double nc=laser0.w0*laser0.w0*mass/4./Pi/q/q;  //Volumetric density
	double n0=density*nc;//15*nc;					//Volumetric density
	
	//Parameter for the current calculation
	double current_factor=nc*thickness/deltag/np; 

	//Limits of the material sheet
    int ileft=int(xgleft/deltag);
	int iright=int(xgright/deltag);
	
	//int na=3;
	
	double qm=-q/mass;
	double ae=0.5*qm*deltat;

		
	particle *elec;
	elec=new particle[np];
	cell1d *cell;
	cell=new cell1d[ndg];
	
	for(int l=0; l<np; l++)
		elec[l].zn=density;
	
	for(int l=0; l<ndg; l++)
	{	
		cell[l].vx.x=(l)*deltag;	
		cell[l].ct.x=(l)*deltag+(deltag/2);
	}
	
	double xg;//,deltag;//deltat,tmax,gamma;
	double Epx,Epy,Epz,Bpx,Bpy,Bpz;
	
	//FFT variables
	double sp1[ndt],sp2[ndt],sp3[ndt],sp4[ndt];//,raux;
	 
	//Particle parameters
	vector<double> qp(np,0.),mp(np,0.);
	
	int iflagion=0;
	int labelx=1, labely=2, labelz=3; //labels used as arguments for the function CalculaJ
	

		
    //Printing parameters in a file
    //fprintf(fo1,"PARAMETROS DEL CAMPO\nE0=%e\nw0=%e\nT0=%e\nlambda=%e\nnw=%d\n\n", laser0.E0, laser0.w0, 2*Pi/laser0.w0, laser0.lambda, laser0.nw);
   // fprintf(fo1,"PARAMETROS DEL MATERIAL\nxgleft=%e\nxfright=%e\nnc=%e\nn0=%e\nileft=%d\niright=%d\n\n", xgleft, xgright, nc, n0, ileft, iright);
   // fprintf(fo1,"\nPARAMETROS DEL CALCULO\nLg=%e\ntmax=%e\nndg=%d\nndt=%d\ndeltag=%e\ndeltat=%e\n\n", lg, tmax, ndg, ndt, deltag, deltat);
   // fprintf(fo1,"PARAMETROS DE LAS PARTICULAS\nnp=%d\nna=%d\nqm=%e\n\n", np, na, qm);
   // fprintf(fo1,"n    t      Eyreflejado    Eytransm    Eyincid    Exleft    Exright    Roleft    Roright\n", np, na, qm);
    //fprintf(fo2,"n(tiempo)   j(num de particula) xp     vpx            vpy         vpz\n");

	/******************************************************/
	/* Initialize variables
	/* Calculate Densities, Fields and position at the beggining
	/******************************************************/
	

	if (iflagion==0)
	{
		npion=np;
		for(int j=0;j<np;j++)
		{
			//qp[j]=-q*n0;
			//mp[j]=mass*n0;
			qp[j]=-q*n0*thickness/np;
			mp[j]=mass*n0*thickness/np;

		}
	}
	else
	{
		npion=0;
	}
	
	//printf("deltag=%f\n", deltag);

	//Initialize particles at rest.	
	if(iflagion==0)
	{
		for(int i=0;i<npion;i++)
		{
			elec[i].xp=thickness/npion*1.*(i+1)+xgleft;
			elec[i].xpold=thickness/npion*1.*(i+1)+xgleft;
			double gamma=1./sqrt(1.-(elec[i].vpx*elec[i].vpx+elec[i].vpy*elec[i].vpy+elec[i].vpz*elec[i].vpz)/lightC_au/lightC_au);
			elec[i].vpx=elec[i].vpx*gamma;
			elec[i].vpy=elec[i].vpy*gamma;
			elec[i].vpz=elec[i].vpz*gamma;
		}
	}//end if
	
	//Calculate the velocities;
 	if(ChargeCalculation==0)
		CalculaRo(elec,cell,qp, deltag, na, np, ndg, fp0);
	else
		for(int k=0; k<np; k++)
			Charge_deposition (cell, elec[k], deltag, current_factor);
	
	if(CurrentCalculation==0)
	{
		CalculaJ(labelx,elec,cell,qp, deltat, deltag, na, np, ndg, fp0, fp0); //labelx is 1 and means 'x'
		CalculaJ(labely,elec,cell,qp, deltat, deltag, na, np, ndg, fp0, fp0); //labely is 2 and means 'y'
		CalculaJ(labelz,elec,cell,qp, deltat, deltag, na, np, ndg, fp0, fp0); //labelz is 3 and means 'z'
	}
	else
	{	printf("Calculo como Lichters la corriente. npion=%d\n", npion);
		for(int k=0; k<npion; k++)
		{	//printf("Llamada a la funcion corriente\n");
			Current_deposition (cell, elec[k], deltag, current_factor);
		}
	}		
			
	
	
	//for(int k=0; k<ndg; k++)
	//	fprintf(fo2,"%d\t%e\n", k, cell[k].vx.Jx);
	
	//Density of positive charge
	if(iflagion==0)
	{
		for(int i=0;i<ndg;i++)
			cell[i].ct.Ropx=-cell[i].ct.Rox;
	}
	else
	{
		for(int i=0;i<ndg;i++)
		{
			cell[i].ct.Ropx=0.;
		}		
	}//end if
	

	
	/******************************************************/
	/* Central part of the program
	/* Time and particle loop 
	/******************************************************/

	for(int n=0;n<ndt;n++)
	{
		if(n%100==0) 
		{
			printf("n=%d\n", n);
			//if(n==100)
			fprintf(fo2,"n=%d\t", n);
			fprintf(fo3,"n=%d\t", n);
			for(int k=0; k<ndg; k++)
			{	fprintf(fo2,"%e\t", cell[k].vx.Jx);
				fprintf(fo3,"%e\t", cell[k].ct.Jy);
			}
			fprintf(fo2,"\n");
			fprintf(fo3,"\n");
		}

		
		double t=n*deltat;
		for(int j=0;j<npion;j++)
		{
			
            interpola(elec[j].xp,cell,Epx,Epy,Epz,Bpx,Bpy,Bpz,t,ndg,deltag,na, fp0);
			
			//Boris algorithm to solve the equation of motion
            ParticleMover(elec[j],Epx,Epy,Epz,Bpx,Bpy,Bpz,deltat,ae);
			
			//Printing positions and velocities
			//if(j%40==0 && n%(ndt/1024)==0)
			  //        fprintf(fo2,"%d\t%d\t%e\t%e\t%e\t%e\n",n, j, xp[j], vpx[j], vpy[j], vpz[j]);
		
		}//End of the Boris algorithm to solve the equation of motion
		if(ChargeCalculation==0)
		{
			CalculaRo(elec,cell,qp, deltag, na, np, ndg, fp0);
		}		
		else
		{	
			for (int l=0;l<ndg;l++)
				cell[l].ct.Rox=0.;
	
			for(int k=0; k<np; k++)
				Charge_deposition (cell, elec[k], deltag, current_factor);
		}
		if(CurrentCalculation==0)
		{
			CalculaJ(labelx,elec,cell,qp, deltat, deltag, na, np, ndg, fp0, fp0); //labelx is 1 and means 'x'
			CalculaJ(labely,elec,cell,qp, deltat, deltag, na, np, ndg, fp0, fp0); //labely is 2 and means 'y'
			CalculaJ(labelz,elec,cell,qp, deltat, deltag, na, np, ndg, fp0, fp0); //labelz is 3 and means 'z'
		}
		else
		{
			for (int l=0;l<ndg;l++)
			{
				cell[l].vx.Jx=0.;
				cell[l].ct.Jy=0.;
				cell[l].ct.Jz=0.;
			}
			for(int k=0; k<npion; k++)
			{	//if(n<2) printf("Llamada a la funcion corriente, %d\n", n);
				Current_deposition (cell, elec[k], deltag, current_factor);
			}
			
		}
		//printf("Rox[205]=%f\tRoy[205]=%f\tRoz[205]=%f\tRox[298]=%f\tRoy[298]=%f\tRoz[298]=%f\n", cell[297].ct.Rox, cell[297].ct.Roy, cell[297].ct.Roz, cell[298].ct.Rox, cell[298].ct.Roy, cell[298].ct.Roz);
		//printf("%d, Rox[297]=%e\tJx[297]=%e\tJy[297]=%e\tJz[297]=%e\tRox[298]=%e\tJx[298]=%e\tJy[298]=%e\tJz[298]=%e\n",n, cell[297].ct.Rox, cell[297].vx.Jx, cell[297].ct.Jy, cell[297].ct.Jz,cell[298].ct.Rox, cell[298].vx.Jx, cell[298].ct.Jy, cell[298].ct.Jz);
		if(n==1)
		{	float Q=0.;
			for(int i=0; i<ndg; i++)
			{	//printf("[%d]=%e\n",i, cell[i].ct.Rox);
				Q+=cell[i].ct.Rox;
			}
			printf("QTOTAL=%e\n", Q);
		}
		
		
		 //if(FieldsCalculation==1)
        //{                      
            /******************************************************/
    		/* Calculate Ex, Ey, Ez, Bx, By, Bz, directly
    		/* By solving the Maxwell equations 
    		/******************************************************/
            
            CalculaCamposInternos(cell, deltat, xgleft, laser0, ndg, deltag, fo1);
       	    for(int k=0; k<ndg; k++)
            {   cell[k].vx.Ey+=E0y(k*deltag,t,xgleft,laser0);
				
                cell[k].vx.Bz+=B0z(k*deltag,t,xgleft,laser0);
            }
			//printf("\tE0y[205]=%e\n", E0y(205*deltag,t,xgleft,laser0));
			//printf("Ex[297]=%f\tEy[297]=%f\tEz[297]=%f\tEx[298]=%f\tEy[298]=%f\tEz[298]=%f\n", cell[297].vx.Ex, cell[297].vx.Ey, cell[297].vx.Ez, cell[298].vx.Ex, cell[298].vx.Ey, cell[298].vx.Ez);

        //}
        //else
		//{
		    /******************************************************/
    		/* Calculate Phi, Ax, Ay and Az
    		/* By solving the wave equation 
    		/******************************************************/
    		
           	/*CalculaPotenciales(phi, Ax, Ay, Az, Axold, Ayold, Azold, DAx, DAy, DAz, Jx, Jy, Jz,IntJx, IntJy, IntJz, Ro, IntRo, Rop, Roold, Ropold, IAxM, IAyM, IAzM, IAxP, IAyP, IAzP, IPhiM, IPhiP, deltat, fp0);
    		
            CalculaEx(Ex,phi,DAx,t+deltat,deltat, fp0);
    		CalculaEy(Ey,DAy,t+deltat,deltat,xg,xgleft,laser0);
    		CalculaEz(Ez,DAz,t+deltat,deltat);		
    		CalculaBy(By,Az,t+deltat,deltat);
    		CalculaBz(Bz,Ay,t+deltat,deltat,xg,xgleft,laser0);
        }*/		
	
		

		
		
		if(iflagion!=0)
		{
			//Ionization();
			//Calculap0();
		}
		
		
		//Calculate the transmitted and reflected field.
		sp1[n]=cell[int(((xgleft-laser0.lambda)/deltag)-1)].vx.Ey-E0y( (xgleft-laser0.lambda)-1,t+deltat,xgleft,laser0);//wp,kp,E0,k0,w0);
		sp2[n]=cell[int(((xgright+laser0.lambda)/deltag)-1)].vx.Ey;
		
		//Printing in fo1, output1
        if(n%(ndt/1024)==0)
        {
                   fprintf(fo1, "%d\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n", n, t/2/Pi*laser0.w0, sp1[n], sp2[n],E0y(xgleft-laser0.lambda,t+deltat,xgleft,laser0), cell[ileft].vx.Ex, cell[iright].vx.Ex, cell[ileft].ct.Rox, cell[iright].ct.Rox);
    }
	}//end time loop
	  
	
	fclose(fo1);
    fclose(fo2);
	fclose(fo3);

	
	
}

void f1(grid a)
{
	cout << a.dx1;
}
